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RESUMEN 

Hemos llevado a cabo simulaciones del tipo Automata Celular para un modelo de 
un polielectrolito en dilucion infinita, para reproducir de manera cualitativa sus 
propiedades conformacionales. Nuestros resultados predicen las llamadas estructuras 
de collar de perlas, las cuales se comparan bien con simulaciones de Dinamica Molec- 
ular mas elaboradas y costosas. 



CELLULAR AUTOMATA SIMULATION OF THE SPATIAL CONFORMATIONS 

OF POLYELECTROLYTES. 

ABSTRACT 

We carried out a Cellular Automata simulation of a model polyelectrolyte solution at 
infinite dilution, in order to reproduce qualitatively its conformational properties. Our 
results predict the so called pearl necklace structures, which compare favorably with 
the more elaborated and costly Molecular Dynamics simulations. 

PACS numbers: 



I. INTRODUCTION 

Polyelectrolyte (PE) solutions are systems widely 
studied since they show properties that are of funda- 
mental interest for applications in health science, food 
industry, water treatment, surface coatings, oil industry, 
among other fields. In fact, one of the problems found in 
genetic engineering in the appearance of conformational 
changes of the ADN molecule, which is a charged 
polyelectrolyte. Q . 

Here we study an infinite dilution polyelectrolyte 
solution, so that, the interaction among polyelectrolyte 



macromolecules are negligible. We model the poly- 
electrolyte as having dissociable functional groups that 
give rise to charged sites and counter-ions in aqueous 
solution. The long range interactions arising from these 
multiple charges are responsible for their macroscopic 
complex properties, which can not be explained by 
regular polymer theories. The spatial structures of these 
materials in solution have been studied extensively, 
particularly with a scaling theory[l|, Q that are not 
appropriate for highly charged PE. The first simulations 
carried out for a single chain predicted the formation 
of groups of monomers, as the fraction of charged 
monomers increased. Such structures are known as 
pearl necklaces. The size of such pearls and the distance 
between then is determined by the balance between the 
electrostatic repulsion and steric effects. 
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These pearl necklace structures have also been found 
in Molecular Dynamics (MD) simulations^!,!^.!!!!, 
fTol ] . In this paper we are interested in the application of 
the much simpler Cellular Automata simulation to char- 
acterize the main features of a polyelectrolytc that could 
be responsible for such conformations. The complete sim- 
ulation of this complex system requires the description 
of a model in terms of potential or forces. In the MD 
simulations of Limbach and Holm [i[ , the monomers are 
connected along the chain by the finite extendible non- 
linear elastic (FENE) bond represented by the potential 
energy. 




U E {r) = -±K E I%hx(l-^ 



(1) 



where r is the distance between two bonded monomers, 
Ke, is the elastic bond constant, a is the monomer diam- 
eter, ks is Boltzmann's constant, T is the absolute tem- 
perature and the parameter Rq represents the maximum 
extension of the bond between two neighbor monomers. 
Two charged sites i and j , with charges eqi and eqj , a dis- 
tance Tij apart, interact with the electrostatic Coulomb 
potential 



Uc{r%j) = k B T 



\Bqiqj 



(2) 



This potential is weighed by the Bjerrum length Xb = 
e 2 /[47re e s /csT], where e s and e D are the solvent permi- 
tivity and the vacuum permitivity respectively, e is the 
electric charge unit. The parameter As is a measure of 
the strength of the electrostatic force as compared to the 
kinetic energy. The length ratio u/Xb is a measure of 
the reduced temperature T/[e 2 /(47re e s fcscr)]. 

The short range and van der Waals interaction between 
any two particles or monomers is represented in the MD 
simulation by a typical truncated Lennard- Jones poten- 
tial 




~t~ ^cut Tij ^ Rc 
TiA > R r 



where e is the potential energy well depth and e cut is the 
cut off energy. This potential prevents the superposition 
of the bonding monomers. Counter-ions interact via a 
purely repulsive LJ interaction with R c — 2 1 / 6 er. 



II. CELLULAR AUTOMATA MODEL 

Even though in the Cellular Automata simulation we 
do not use any form of potential energies or forces in 
an explicit manner, the rules for the movement of the 
different particles must be inspired on a model defined 
in terms of such potentials. We therefore establish 
our rules based on the essence of the previous three 



FIG. 1: Celular Automata Model. The dark bonded circles 
represent the charged monomers, the gray circles are the neu- 
tral monomers and the open circles are the counter-ions. The 
shaded area corresponds to the range volume Xb = 3<r 



potentials. 



A. Polyelectrolyte Chain Construction Rules 

The polymer is constructed by placing the monomers 
in a three dimensional cubic network of side L and 
volume L 3 . Each cell then has 26 neighbors and rep- 
resents a monomer with a monovalent charge +1, — 1, 
or, for a neutral monomer, 0, as depicted in Fig.(l) in 
two dimensions. Out of a total of N monomers in the 
chain, it is assumed that a given fraction / is charged. 
The polymer is then constructed by randomly binding 
consecutive sites in the network. Each monomer could 
be charged or uncharged, with a distribution chosen 
randomly. A key step on the construction of the polyion 
is the spatial location of the dissociated counter-ions. 
We place the counter-ions also randomly in free cells 
in the volume around the charged monomer within a 
distance As, that is, in a volume (2Ab) 3 centered on 
the charged site. The use of the Bjerrum parameter, 
which is related to the quality of the solvent, ensures 
the conservation of the total electroneutrality but gives 
a spatial distribution of counter-ions around the charged 
sites. 

So, each monomer i of the system is represented in 
a NxA matrix where each element M(m x ,m y ,m z ,qi) 
indicates the i th polymer with charge that could be 
— 1, +1 or 0, at the positions x, y, z given by the cell label 
[m x , m v , m z ]. The counter-ions with opposite charge are 
represented by a similar f NxA matrix Mc- 

For simplicity we chose monomers with dissociable 
groups that give a site with a positive charge. We then 
set the following displacement rules for the different 
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particles: 



B. Displacement Rules 
Neutral Monomer Particle 

1. Locate the unoccupied nearest neighbor sites. The 
new position where a move could be acceptable 
are those where it does not superimposed with any 
other particle in the system and where no bond is 
broken. 

2. Count the amount of monomer particles around 
the current position and around every unoccupied 
neighbor, within a cube of volume (2L) 3 centered 
on it. 

3. Move the test neutral monomer to the position that 
has the higher amount of monomers around it, in- 
cluding the current one if it were the case. 

I 




FIG. 2: Final poly electrolyte configuration for different charge 
conformations appear as the charge fraction increases. 



III. RESULTS AND DISCUSSION 



Positively Charged Monomer Particle 

1. As before, locate the unoccupied acceptable nearest 
neighbor sites. 

2. Count the amount of charge around the current 
position and around every unoccupied neighbor, 
within a cube of volume (2A_b) 3 centered on it. 

3. Move the test positive charged monomer to the po- 
sition that has the lowest positive charge around it, 
including the current one if it were the case. 

Negatively charged counter-ion 

1. Move randomly to an unoccupied site within a cube 
of volume (2As) 3 centered on the accepted new po- 
sition of the corresponding positive monomer in the 
polymeric chain. 




fractions f: A) 0.1 B) 0.3 C) 0.5 D) 0.7. The pearl necklace 



I 

We denote the position of monomer i with r» and the 
distance between two particles i and j with ry. The 
center of mass for the chain is then R s = jfJ2iLi r i- 
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FIG. 3: Evolution of the Radius of Gyration R g as a function of time steps t, for N = 100, As/cr = 3 and f 
corresponds to different random charged sites distributions, obtained with different seeds. 



0.5. Each curve 



and the center of mass coordinates are a;, = r.; — R s . 
A parameter that is useful in the study of the spatial 
conformations of a polymer is the radius of gyration Rq 1 
defined as 



1 N 

i=l 



(4) 



According with our construction and movement rules, 
we can vary the length of the chain L, or the number 
of monomers N, satisfying L — Na and the number 
of charged monomers, N c — JN. We also take as an 
independent variable the parameter Xg, which deter- 
mines the number of cells, of size a, where the range of 
the electrostatic attraction between a charged monomer 
and contra-ion extends. The charge distribution of the 
sequence of charged and uncharged monomers is de- 
termined randomly depending on the initial random seed. 

In Fig. (2) we show some of the equilibrium confor- 
mations obtained for a polyelectrolyte with N = 100 
monomers, with As = 3<r, for several charge fractions. 
The line represents the polyelectrolyte, the filled black 
bonded circles represent the charged monomers. The 
counter-ions in solution are represented by open non- 
bonded circles. For clarity, the neutral monomers are 



not shown. For a low charge fraction of / = 0.1, the 
polyion presents an elongated string appearance. As 
the charge fraction is increased the polyion contracts 
and some groups of monomers tend to form clusters, so 
that, already for / = 0.5, it shows the locally collapsed 
structures known as pearl necklace. These results 
are very similar to those obtained by the Molecular 
Dynamics simulations of Limbach and Holm 0, 0] 



It is important to notice that for a given total charge, 
determined by the fraction /, different distributions of 
the charged monomers give different conformations. To 
study this behavior, we have carried out simulations 
with several initial seeds for N = 100 monomers, with 
a fixed fraction / = 0.5 and a range parameter of 
Ab/<t = 3. In Fig. (3), we show the temporal evolution of 
the radius of gyration R g . In all cases the conformations 
change from the initial given R g value to a plateau 
value that correspond to the equilibrium structures. 
Fig. (3) clearly show that the plateau R g values, and 
thence, the final conformations depend on the charge 
distribution. In Fig. (4) we show snapshots of final 
structures corresponding with the different seeds of 
Fig.(3). 



J 




As we can see from Fig. (4), the formation of the pearl 
necklace structures is independent of the charge sites dis- 
tribution, for a given / and Xb- These clusters seem to 
be stabilized by the counter-ions as a consequence of the 
electroneutrality condition that we force to be satisfied. 
This is so because in our model the degree of freedom of 



the mobile counter-ions is much higher that that of the 
monomers tied to the chain. The strong repulsions that 
originate by the formation of clusters of neutral and posi- 
tively charged monomers is compensated by the counter- 
ion cloud that forms around it. 



In order to study the reproducibility of the configura- 
tions found, we carried out a large number of simulation 
runs, for fixed values of TV = 150, / = 0.5, A^/cr = 3 
and 9 and for the same initial charge distribution. In 
Fig. (5) we show the histograms for the frequency with 



which a given value of R g appears. We can see that the 
distribution of the R g is very closely a gaussian with a 
reasonably low dispersion of less than a 10% about the 
mean value. 



In Fig.(5B) we can see that similar distribution is ob- 



tained when the range parameter Ab/ct increases to a 
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value of 9. We have further tested the effect of the pa- 
rameter Xb by generating structures for different values 
of it. In Fig. (6) we show some equilibrium conformations 
for \b/o equal to A) 3,B) 6,C) 9 and D) 12. Here we use 
a large charge fraction / = 0.5 and N = 100 and we used 
the same initial charge distributions in all cases. As we 
can observe as Bjerrum length increases the final poly- 
electrolyte structures become more compact. The num- 
ber or pearls or conglomerates is higher for the lower 
values of As, a result similar to that obtained by MD 
simulations of Limbach y Holm pl|-[loj 

IV. CONCLUSIONS 

With the simple technique described here, we were able 
to reproduce the complex structure of model polyelec- 
trolytes that fare very well with those predicted by the 
more sophisticated Molecular Dynamic and Monte Carlo 
simulations 0, 0| . We even predict situations with single 



conglomerates and with pearl necklace type conglomer- 
ates. We thus show the potentiality of the Celular Au- 
tomata in the simulation of the trends in the formation 
of the various types of spatial conformations of polyclcc- 
trolytes. We remark on the importance of the charge 
distribution once the fractional charge is fixed. 
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FIG. 6: Polyelectrolyte conformations for several values of the reduced Bjerrum parameter \b/o~ A) 3,B) 6,C) 9 and D) 12. 
Here f = 0.9 and N = 100 and we used the same initial charge distributions 



